#!/usr/bin/python

from numpy import *
from math import *
from cmath import exp
import Gnuplot
from scipy.interpolate import UnivariateSpline

from band_struct import *

def Linear_Sigma( ens, nb ):
#
# Parameters a, b
  a = -0.0
  b = -0.0
  c = -0.0
#  a = -0.0
#  b = -0.0
# Parameters Re1, Re2
  Re1 =-0.0
  Re2 = 0.0
  Re3 =-0.0

  assert nb == 4, "  nb must be equal 4"

  ne = len(ens)

  Sigma = zeros( (ne, nb, nb), dtype = complex )

  Sig3_cl = zeros( (3, 3), dtype = complex )
  Sig3_site = zeros( (3, 3), dtype = complex )

# Sigma in the cluster basis
  Sig3_cl[0, 0] = a
  Sig3_cl[1, 1] = a
  Sig3_cl[2, 2] = b

# Transform to site basis
  hloc = array( [ [ 0, 1, 1 ], \
                  [ 1, 0, 1 ], \
                  [ 1, 1, 0 ] ], dtype = float )

  [eigv, eigvecs] = linalg.eigh(hloc)
  
  Sig3_site = dot( eigvecs, dot( Sig3_cl, eigvecs.T.conj() ) )

# Sigma in the cluster basis
  for ie in range(ne):
    Sigma[ie, 0:3, 0:3] = Sig3_site * ens[ie]
    Sigma[ie, 3, 3] = c * ens[ie]

# Sigma in the cluster basis
  Sig3_cl[0, 0] = Re1
  Sig3_cl[1, 1] = Re1
  Sig3_cl[2, 2] = Re2

# Additional constant energy shift (HF part)
  Sig3_site = dot( eigvecs, dot( Sig3_cl, eigvecs.T.conj() ) )

  for ie in range(ne):
    Sigma[ie, 0:3, 0:3] += Sig3_site
    Sigma[ie, 3, 3] += Re3

# 
  return Sigma

##################################################################
#
# Sigma_of_k
#
##################################################################
def Sigma_of_k_kag_3x3(Efermi, nb, sigfname = None):
  """
  Prepares and returns a function evaluating the k-dependent
  self-energy.
  """

  Sigma_fun = Read_Sigma(sigfname, nb)

  print " Im Sigma( 0.0 ):"
  print imag( Sigma_fun(0.0) )

  def Sigma_k_fun( en, kv ):
#    tmp = zeros( en.shape + (3, 3), dtype = complex )
#    tmp[:, 0, 0] = 0.01j - Efermi*0
#    tmp[:, 1, 1] = 0.01j - Efermi*0
#    tmp[:, 2, 2] = 0.01j - Efermi*0
#    return tmp
#
# Periodized Sigma
#
    k1 = kv[0]
    k2 = kv[1]
    k21 = k2 - k1

    Sigma = Sigma_fun(en)
    Sigma_k = array( Sigma )

# Self-energy periodization
    Sig_0 = ( Sigma[:, 0, 0] + Sigma[:, 1, 1] + Sigma[:, 2, 2] ) / 3 - Efermi
    Sig_1 = ( Sigma[:, 0, 1] + Sigma[:, 0, 2] + Sigma[:, 1, 2] ) / 3

    Sigma_k[:, 0, 0] = Sig_0
    Sigma_k[:, 1, 1] = Sig_0
    Sigma_k[:, 2, 2] = Sig_0
    Sigma_k[:, 1, 0] = 0.5 * Sig_1 * (1.0 + exp( -2.0j * k1 ) )
    Sigma_k[:, 2, 0] = 0.5 * Sig_1 * (1.0 + exp( -2.0j * k2 ) )
    Sigma_k[:, 2, 1] = 0.5 * Sig_1 * (1.0 + exp( -2.0j * k21 ) )
    Sigma_k[:, 0, 1] = 0.5 * Sig_1 * (1.0 + exp(  2.0j * k1 ) )
    Sigma_k[:, 0, 2] = 0.5 * Sig_1 * (1.0 + exp(  2.0j * k2 ) )
    Sigma_k[:, 1, 2] = 0.5 * Sig_1 * (1.0 + exp(  2.0j * k21 ) )


# Cumulant periodization
#    M = zeros( Sigma.shape, dtype = complex )
#    M_k = zeros( Sigma.shape, dtype = complex )
#    umat = eye( nb, dtype = complex )
#
#    for ie, e in enumerate( en ):
#      M[ie, : :] = linalg.inv( (e + Efermi) * umat  - Sigma[ie, :, :] )
#
#    M_0 = ( M[:, 0, 0] + M[:, 1, 1] + M[:, 2, 2] ) / 3
#    M_1 = ( M[:, 0, 1] + M[:, 0, 2] + M[:, 1, 2] ) / 3
#
#    M_k[:, 0, 0] = M_0
#    M_k[:, 1, 1] = M_0
#    M_k[:, 2, 2] = M_0
#    M_k[:, 1, 0] = 0.5 * M_1 * (1.0 + exp( -2.0j * k1 ) )
#    M_k[:, 2, 0] = 0.5 * M_1 * (1.0 + exp( -2.0j * k2 ) )
#    M_k[:, 2, 1] = 0.5 * M_1 * (1.0 + exp( -2.0j * k21 ) )
#    M_k[:, 0, 1] = 0.5 * M_1 * (1.0 + exp(  2.0j * k1 ) )
#    M_k[:, 0, 2] = 0.5 * M_1 * (1.0 + exp(  2.0j * k2 ) )
#    M_k[:, 1, 2] = 0.5 * M_1 * (1.0 + exp(  2.0j * k21 ) )
#
#    for ie, e in enumerate( en ):
## N.B. Efermi is subtracted!
#      Sigma_k[ie, : :] = e - linalg.inv( M_k[ie, :, :] )

    return Sigma_k

  return Sigma_k_fun


def StoreData(base):
#
# Store data in 3 files:
#   <base>.eig - eigenvalues
#   <base>.plo - projectors
#   <base>.kpoints - k-points
#
#
# Filenames
#
  def transf_apply( A, U ):
    return dot( U.conj().T, dot( A, U ) )

  FeigName = base + ".eig"
  FmatelName = base + ".matel"
  FploName = base + ".plo"
  FkptsName = base + ".kpoints"
  FklatName = base + ".klat"

  nk = 41
  nkz = 1
  ne = 501
  nb = 3
  no = 3

  thopz = 0.0

  dk = pi / ( nk  )
  kr1 = linspace(-pi/2, pi/2 - dk, nk)
  kr2 = linspace(-pi/2, pi/2 - dk, nk)
# A k-mesh extended to BZ boundaries
#  kr1 = linspace(-pi/2, pi/2, nk)
#  kr2 = linspace(-pi/2, pi/2, nk)

  nktot = nk*nk*nkz

  kweight = 1.0/nktot

  b1 = array([1, -1/sqrt(3)])
  b2 = array([0,  2/sqrt(3)])

  eigk = zeros( (nk, nk, nb), dtype = float )
  hamk = zeros( (nk, nk, nb, nb), dtype = complex )
  Uk = zeros( (nk, nk, nb, nb), dtype = complex )
  matel = zeros( (2, nb, nb), dtype = complex )
#
# Obtain e(k) for each k-point and store to the .eig-file
#
  feig = open(FeigName, "wt")
  fmatel = open(FmatelName, "wt")
  fplo = open(FploName, "wt")
  fkpts = open(FkptsName, "wt")
  fklat = open(FklatName, "wt")

  fplo.write("%10d  %8d  %8d  %8d  %8d\n\n\n"%(nb, no, nktot, 1, 1))

  fklat.write("%8d  %8d  %8d\n"%(nk, nk, nkz))
  fklat.write("%15.10f  %15.10f  %15.10f\n"%(b1[0], b1[1], 0.0))
  fklat.write("%15.10f  %15.10f  %15.10f\n"%(b2[0], b2[1], 0.0))
  fklat.write("%15.10f  %15.10f  %15.10f\n\n"%(0.0, 0.0, 1.0))

  ikp = 0

  for ik2 in range(nk):
    for ik1 in range(nk):
 
      k1 = kr1[ik1]
      k2 = kr2[ik2]

      [ek, eigvecs] = kag_eig(k1,k2)
#      [ek, eigvecs] = tri4_eig(k1,k2)
#      [ek, eigvecs] = tri_eig(k1,k2)

      Uk[ik1, ik2, :, :] = eigvecs

      kv = k1*b1 + k2*b2
#
# z-component
#
      for ikz in range(nkz):
        ikp += 1

        if nkz > 1:
          kz = float(ikz - nkz/2) / ( nkz / 2 ) * pi
          if ik1 + ik2 == 0: print "Kz:", kz
        else:
          kz = 0.0

        ekz = -2 * thopz * cos(kz)

        fkpts.write("%20.15f %20.15f %20.15f %20.15f\n"%(kv[0], kv[1], kz, kweight))
        fklat.write("%20.15f %20.15f %20.15f\n"%(kv[0], kv[1], kz))


#        hamk[ik1, ik2, :, :] = transf_apply( tri4_ham(k1, k2), eigvecs )
#        hamk[ik1, ik2, :, :] = tri4_ham(k1, k2)
        hamk[ik1, ik2, :, :] = kag_ham(k1, k2)

        for ib in range(nb):
          feig.write("%10d %9d %29.20f\n"%(ib+1, ikp, ek[ib] + ekz))

        fplo.write("\n")

        for io in range(no):
          for ib in range(nb):
            proj = eigvecs[io, ib]
            fplo.write("%30.20f %30.20f\n"%(real(proj), imag(proj)))
  
  feig.write("\n")

  kbas = vstack( (b1, b2) )

  print "  Basis in recirpocal space:", kbas
  print "  Basis in direct space:", linalg.inv(kbas)

#  velk = Numerical_Velocities( eigk, kbas )
# Atom positions in the cluster
  sy = sqrt(3.0) / 2
  svec = array( [[ 0.0, 1.0, 0.5, 1.5 ], \
                 [ 0.0, 0.0,  sy,  sy ]] ) 

#  print "!!!!!!!!!!!!!!!!!"
#  print " svec is set to zero!"
#  print "!!!!!!!!!!!!!!!!!"
#
#  svec[:] = 0.0
  melk = Numerical_Matrix_Elements( hamk, svec, kbas, [dk, dk] )

#
# Matrix elements
#
  print
  print "  Storing matrix elements..."
  ikp = 0
  averr = 0.0
  for ik2 in range(nk):
    for ik1 in range(nk):
      k1 = kr1[ik1]
      k2 = kr2[ik2]

#      [vx, vy] = tri_vel( k1, k2 )

#      err = (vx - velk[0, ik1, ik2, 0])**2 + (vy - velk[1, ik1, ik2, 0])**2
#      averr += err
#      assert err < 1.0e-3/nk, \
#             "  Something is wrong: %d  %d"%(ik1, ik2)
      
      ikp += 1

      for i in range(2):
        matel[i, :, :] = transf_apply( melk[i, ik1, ik2, :, :], Uk[ik1, ik2, :, :] )
#        matel[i, :, :] = melk[i, ik1, ik2, :, :]

      for ib2 in range(nb):
        for ib1 in range(nb):
          for i in range(3):
            if i < 2:
              vr = real(matel[i, ib1, ib2])
              vi = imag(matel[i, ib1, ib2])
            else:
              vr = 0.0
              vi = 0.0

            fmatel.write("%10d %9d %9d %9d %29.20f %29.20f\n"\
                         %(ib1+1, ib2+1, ikp, i+1, vr, vi))

      fmatel.write("\n")

#  print "  Average velocity error:", averr / nk / nk

  feig.close()
  fmatel.close()
  fkpts.close()
  fklat.close()
  fplo.close()

#
# Evaluates Fermi velocities numerically
# as the gradient of eigenvalues
#
def Numerical_Velocities( eigk, kbas ):

  coefs = array( [1, -8, 8, -1], dtype = float) / 12.0

  abas = linalg.inv(kbas)

  nkx = eigk.shape[0]
  nky = eigk.shape[1]
  nb = eigk.shape[2]

  A = array( [abas[:,0] * nkx, abas[0:,1] * nky] ).T / 2 /pi

  velk = zeros( (2, nkx, nky, nb), dtype = float )

  p1 = zeros( 4, dtype = int )
  p2 = zeros( 4, dtype = int )
  for ik2 in range(nky):
    for ik1 in range(nkx):
      p1[0] = (ik1 - 2 + nkx) % nkx
      p1[1] = (ik1 - 1 + nkx) % nkx
      p1[2] = (ik1 + 1) % nkx
      p1[3] = (ik1 + 2) % nkx

      p2[0] = (ik2 - 2 + nky) % nky
      p2[1] = (ik2 - 1 + nky) % nky
      p2[2] = (ik2 + 1) % nky
      p2[3] = (ik2 + 2) % nky

      for ib in range(nb):
        v1 = sum( eigk[p1, ik2, ib] * coefs )

        v2 = sum( eigk[ik1, p2, ib] * coefs )

        velk[0, ik1, ik2, ib] = A[0,0] * v1 + A[0,1] * v2
        velk[1, ik1, ik2, ib] = A[1,0] * v1 + A[1,1] * v2
  
  return velk

#
# Evaluates Fermi velocities numerically
# as the gradient of eigenvalues
#
def Numerical_Matrix_Elements( hamk, svec, kbas, dk ):

  coefs = array( [1, -8, 8, -1], dtype = float) / 12.0

  abas = linalg.inv(kbas)

  nkx = hamk.shape[0]
  nky = hamk.shape[1]
  nb = hamk.shape[2]

  A = array( [abas[:,0] * nkx, abas[0:,1] * nky] ).T / 2 /pi

  velk = zeros( (2, nkx, nky, nb, nb), dtype = complex )
  sij = zeros( (2, nb, nb), dtype = float)

  for ib2 in range(nb):
    for ib1 in range(nb):
      sij[:, ib1, ib2] = svec[:, ib1] - svec[:, ib2]

  p1 = zeros( 4, dtype = int )
  p2 = zeros( 4, dtype = int )

  kscal_1 = 2 * pi / nkx / dk[0]
  kscal_2 = 2 * pi / nky / dk[1]

  for ik2 in range(nky):
    for ik1 in range(nkx):
      p1[0] = (ik1 - 2 + nkx) % nkx
      p1[1] = (ik1 - 1 + nkx) % nkx
      p1[2] = (ik1 + 1) % nkx
      p1[3] = (ik1 + 2) % nkx

      p2[0] = (ik2 - 2 + nky) % nky
      p2[1] = (ik2 - 1 + nky) % nky
      p2[2] = (ik2 + 1) % nky
      p2[3] = (ik2 + 2) % nky

      for ib2 in range(nb):
        for ib1 in range(nb):
          v1 = sum( hamk[p1, ik2, ib1, ib2] * coefs ) * kscal_1

          v2 = sum( hamk[ik1, p2, ib1, ib2] * coefs ) * kscal_2

          velk[0, ik1, ik2, ib1, ib2] = A[0,0] * v1 + A[0,1] * v2
          velk[1, ik1, ik2, ib1, ib2] = A[1,0] * v1 + A[1,1] * v2

      velk[0, ik1, ik2, :, :] += 1.0j * sij[0, :, :] * hamk[ik1, ik2, :, :]
      velk[1, ik1, ik2, :, :] += 1.0j * sij[1, :, :] * hamk[ik1, ik2, :, :]
  
  return velk

#
# Test of the matrix elements for the 4-site triangular lattice
#
def Test_tri4_Vel():

  def transf_apply( A, U ):
    return dot( U.conj().T, dot( A, U ) )

  nk = 20
  nb = 3
  no = 1

  dk = pi / ( nk  )
  kr1 = linspace(-pi / 2, pi / 2 - dk, nk)
  kr2 = linspace(-pi / 2, pi / 2 - dk, nk)

  nktot = nk*nk

  kweight = 1.0/nktot

  b1 = array([1, -1/sqrt(3)])
  b2 = array([0,  2/sqrt(3)])

  hamk = zeros( (nk, nk, 4, 4), dtype = complex )
  velnu = zeros( (2, 4, 4), dtype = complex )
  
  for ik2 in range(nk):
    for ik1 in range(nk):
      k1 = kr1[ik1]
      k2 = kr2[ik2]

      hamk[ik1, ik2, :, :] = tri4_ham( k1, k2 )

  kbas = vstack( (b1, b2) )

  sy = sqrt(3) / 2

# Atom positions in the cluster
  svec = array( [[ 0.0, 1.0, 0.5, 1.5 ], \
                 [ 0.0, 0.0,  sy,  sy ]] ) 

  dk_diff = array( [ kr1[1] - kr1[0], kr2[1] - kr2[0] ] )

  melk = Numerical_Matrix_Elements( hamk, svec, kbas, dk_diff )
  
  melk[ abs(melk) < 1e-9 ] = 0.0

  averr = 0.0
  for ik2 in range(nk):
    for ik1 in range(nk):
 
      k1 = kr1[ik1]
      k2 = kr2[ik2]

      melij = tri4_vel( k1, k2 )

      melij[ abs(melij) < 1e-9 ] = 0.0

#      print
#      print " Analytic:"
#      for ib2 in range(nb):
#        for ib1 in range(nb):
#          print ib1, ib2, melij[:, ib1, ib2]
#
#      print
#      print " Numerical:"
#      for ib2 in range(nb):
#        for ib1 in range(nb):
#          print ib1, ib2, melk[:, ik1, ik2, ib1, ib2]
#
#      raw_input("  Press")
      err = 0.0
      for ib2 in range(nb):
        for ib1 in range(nb):
          err += sum( abs( melij[:, ib1, ib2] - melk[:, ik1, ik2, ib1, ib2] )**2 )

      averr += err
      assert err < 1e-5, "  Error is too large: %d, %d, %15.9f"%(ik1, ik2, err)

  print "  Average error:", averr / nk**2

#
# Test of the matrix elements for the Kagome lattice
#
def Test_Kag_Vel():

  def transf_apply( Uconj, A, U ):
    return dot( Uconj, dot( A, U ) )

  nk = 20
  nb = 3
  no = 1

  dk = pi / ( nk  )
  kr1 = linspace(-pi / 2, pi / 2 - dk, nk)
  kr2 = linspace(-pi / 2, pi / 2 - dk, nk)

  nktot = nk*nk

  kweight = 1.0/nktot

  b1 = array([1, -1/sqrt(3)])
  b2 = array([0,  2/sqrt(3)])

  eigk = zeros( (nk, nk, 3), dtype = float )
  velnu = zeros( (2, 3, 3), dtype = complex )
  
  for ik2 in range(nk):
    for ik1 in range(nk):
      k1 = kr1[ik1]
      k2 = kr2[ik2]

      [ek, eigvecs] = kag_eig(k1,k2)
      eigk[ik1, ik2, :] = ek[:]

  kbas = vstack( (b1, b2) )

  velk = Numerical_Velocities( eigk, kbas )
  
  for ik2 in range(nk):
    for ik1 in range(nk):
 
      k1 = kr1[ik1]
      k2 = kr2[ik2]

      [ek, eigvecs] = kag_eig(k1,k2)

      velij = kag_vel( k1, k2 )

# Transform to eigenbasis
      velnu[0, :, :] = transf_apply( eigvecs.conj().T, velij[0, :, :], eigvecs )
      velnu[1, :, :] = transf_apply( eigvecs.conj().T, velij[1, :, :], eigvecs )

      velnu[ abs(velnu) < 1e-14 ] = 0.0

      print " k-point:", ik1, ik2
      print velnu[0, :, :]
      print velnu[1, :, :]

      print 
      print " From band structure:"
      print velk[:, ik1, ik2, :]

      print
      print " Transformed to ij:"
      print transf_apply( eigvecs, diag( velk[0, ik1, ik2, :], 0 ), eigvecs.conj().T )
      print transf_apply( eigvecs, diag( velk[1, ik1, ik2, :], 0 ), eigvecs.conj().T )

      print
      print " Matrix elements ij:"
      print velij[0, :, :]
      print velij[1, :, :]
      raw_input( " Press ")

      kv = k1*b1 + k2*b2
#
def ReadSigmaOld(filename, no):
#
# Reads Sigma defined on real axis from file <filename>
#
  fsig = open(filename)

  line = fsig.readline()
  while not line.strip():
    line = fsig.readline()

  oms = []
  Sig = [[] for i in range(no) ]

  for i in range(no):
    while line.strip():
      sline = line.split()
      if i==0:
        oms.append(float(sline[0]))
      Sig[i].append(complex(float(sline[1]), float(sline[2])))
      line = fsig.readline()

    line = fsig.readline()

  fsig.close()

  return [ array(oms), array(Sig) ]

def Read_Sigma( filename, nb, Inds = ['Do'] ):
#
# Reads Sigma in site representation from the files
# stored by pyDMFT
#
  alloc = False

  SigRe_fun = [ [] for ib in range( nb ) ]
  SigIm_fun = [ [] for ib in range( nb ) ]

  for ind in Inds:
    for i in range(1, nb + 1):
      for j in range(1, nb + 1):
        fileij = filename + "_" + ind + "/" + "%d_%d.dat"%(i, j)

        fsig = open(fileij)

        line = fsig.readline()
        while not line.strip():
          line = fsig.readline()

        oms = []
        Sig = []

        while line.strip():
          sline = line.split()
          oms.append(float(sline[0]))
          Sig.append(complex(float(sline[1]), float(sline[2])))
          line = fsig.readline()

        fsig.close()

        ens = array( oms )
        Sigma = array( Sig )
  
        ib1 = i - 1
        ib2 = j - 1
        SigRe_fun[ib1].append( UnivariateSpline(ens, real(Sigma[:]), k = 3, s = 0) )
        SigIm_fun[ib1].append( UnivariateSpline(ens, imag(Sigma[:]), k = 3, s = 0) )

  def Sigma_fun(en):
    if isscalar(en): en = array( [en] )

    SigwRe = zeros( en.shape + (nb, nb), dtype = float )
    SigwIm = zeros( en.shape + (nb, nb), dtype = float )

    for ib2 in range(nb):
      for ib1 in range(nb):
        SigwRe[:, ib1, ib2] = SigRe_fun[ib1][ib2]( en )
        SigwIm[:, ib1, ib2] = SigIm_fun[ib1][ib2]( en )

    return SigwRe + 1.0j * SigwIm

  return Sigma_fun

##################################################
#
# Beginning of the main program
#
##################################################
graph = Gnuplot.Gnuplot()
V4 =  -0000.2

nk_stacks = 21
ne_stacks = 201

nk_map = 101
ne_map = 101

emin = -1.5
emax =  0.5
en_stacks = linspace(emin, emax, ne_stacks)
en_map = linspace(emin, emax, ne_map)

#mod = BS_Square_L2()
#mod.Generate_Save("square_L2_test", 20)

bs = Band_Structure()

BZ_bas = 2 * pi * array( [[1, 0], \
                          [0, 1]] )

bs.Read_External("square_L2_test", 20, BZ_bas)

assert 0

plotter = Plotter()

#mod = BS_Kag_L3()
mod = BS_Tri4_L3(0.0)

plotter.GraphOutput = "tri4_L3.png"
plotter.GraphTitle = None
plotter.Filename = "test_kpoint.dat"
plotter.onScreen = False

mod.SF_band_structure(0.4, nk_map, en_map, None, plotter, orbs=[0], Stacks3D=False)
#mod.SF_kpoint(0.4, [pi, pi], en_stacks, None, plotter, orbs=[0])
en_fermi = linspace(-0.01, 0.01, 5)
kmax = 6.0 * pi / 3.0
#mod.SF_Fermi_surface(0.4, 51, kmax, en_fermi, None, plotter, unfold=False, orbs=None)

sigfname = "kag0.67/Beta200_V0.10/SigmaRe_site"
#Crystal_Field()
#Density_of_States(6.0196)
#Fermi_Surface(0.580)
#Band_Structure(0.333, True)  # D = 0.56  (U = 0.0)
#Band_Structure(0.580)  # V = -0.0  (x = 0.67)
#Band_Structure(0.500)  # V = -0.5  (x = 0.67)
#Band_Structure(0.357)  # V = -100.0  (x = 0.67)
#Band_Structure(5.14021985704) # D = 0.56  Single-site
#Band_Structure(5.12099077911) # D = 0.56  V = 0.00
#Band_Structure(0.489394674985, True) # D = 0.56  V = 0.00, L3
#Band_Structure(0.38084628118, False) # D = 0.56  V = 0.00, L3, modified CF
#Band_Structure(5.42000280663)  # D = 0.56  V = 0.10
#Band_Structure(6.01962692593) # D = 0.56  V = 0.30
#Band_Structure(0.471532855828, False) # D = 0.56  V = 0.30, L3
#Band_Structure(6.61486)
#Band_Structure(5.18136379555)  # D = 0.62  V = 0.00
#Band_Structure(5.89833981884)  # D = 0.62  V = 0.30, Beta = 100

#Band_Structure(5.23908978742)  # D = 0.67  V = 0.00
#Band_Structure(5.41555134613)  # D = 0.67  V = 0.05
#Band_Structure(5.5929766719)  # D = 0.67  V = 0.10
#Band_Structure(5.61035688621)  # D = 0.67  V = 0.10, Beta = 200
#Band_Structure(5.7669023627)  # D = 0.67  V = 0.15
#KPoint(5.61035688621, [0.0, 0.0], "Ak_G_D0.67_b200_V0.10.sfun") # Gamma point
#KPoint(5.42000280663, [0.0, 0.0], "Ak_G_D0.56_b40_V0.10.sfun") # Gamma point
#KPoint(6.01962692593, [4*pi/3, 2*pi/3], "Ak_K_D0.56_b40_V0.10.sfun") # point K
#KPoint(6.01962692593, [pi, pi], "Ak_M_D0.56_b40_V0.10.sfun") # point M
#KPoint(6.01962692593, [4*pi/3/10.0*9, 2*pi/3/10.0*9], "Ak_GK_9_D0.56_b40_V0.30.sfun") # point at GK
#KPoint(6.01962692593, [pi/10.0*9, pi/10.0*9], "Ak_GM_9_D0.56_b40_V0.30.sfun") # point at GM
#StoreData("kag_t3")
#StoreData("tri4_t3_v_0.0_80")
#Test_tri4_Vel()
#mod = BS_Square_L2()
#mod.Generate_Save("sq_L2s_k10_td0.2_to0.05_de0.0", 10)
#mod = BS_Kag_L3()
#mod.Generate_Save("kag_L3_p1_k20_m", 21)
#mod = BS_Tri4(0.0)
#mod.Generate_Save("tri4_v0.0_k20", 20)
#mod.HMLT_Save("kag_L3_p1_CF150.0_k40", 41)
#mod.HMLT_Save("sq_test", 10)

#print " Reading Sigma..."
#Sigma_fun = Read_Sigma("kag0.56/V0.00/SigmaRe_site", 3)
#print Sigma_fun(array([0,1,2,3,4])).shape

raw_input('Press Enter...')


